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1. INTRODUCTION 

It is well known that boundary element techniques require the 
integration of a characteristic solution over elements which 
define the boundary of the domain. The characteristic solution 
is usually a function of the distance from a "root” node to the 
surface of the element over which the integral is being computed. 
When the "root* node is outside of the element, the integrals may 
be performed in a straight-forward fashion using standard numeri- 
cal integration techniques like Gaussian quadrature. However, 
when the node is in the element, the integral contains an integr- 
able singularity which will not integrate accurately using Gaus- 
sian quadrature. This report will discuss an alternative method 
for performing this integral. 

The method proposed separates the integral of the charac- 
teristic solution into a singular and non-singular part. The 
singular portion is integrated with a combination of analytic and 
numerical techniques while the non-singular portion is integrated 
with standard Gaussian quadrature. The method may be generalized 
to many types of sub-parametric elements. 

This report will consider only the integrals over elements 
containing the root node, and will deal only with the charac- 
teristic solution for linear acoustic problems. However, the 
method described may be generaxized to most characteristic solu- 


tions. 
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2. BEM INTEGRALS OVER AN ELEMENT CONTAINING THE ROOT HOPE 

Suppose the root node is node n (1 < n < m) . To denote dis- 
tance from node n, write r as r . 

n 

e' 31 “ n 

Nj — dA i ■ 1.2, ... r m (1) 

e n 

where 

I A - contribution (i.e., integral) to coefficient of node i. 

- shape function for node i 
A tt - area of element e (region of integration) 
r n - distance from root node to any point in the element 
k - constant from Helmholtz equation 
m - number of nodes in element e 
Several significant properties of the shape functions are: 

1) All {N^i * n) —* 0 as r n — 0 

2) N — 1 as r_ — * 0 

' n n 

m 

3) N n ‘ 1 - L V 

n i-1 1 

i*n 


From these observations we note that the integrand of 1^ is 
singular only when i»n. For all other integrals, — 0 as 
r n — * 0 quickly enough to avoid singularity. Thus, I ir for i*n 
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may be integrated numerically without undo difficulty. 


Por i»n an integrable singularity exists. First, rewrite 

the integral I as: 
n 


*n - /* N n dA 

e n 


-jkr 


- J 


n 


A e r n 


dA - 


ro 

t I 
i-1 
i*n 


CD 


We may therefore break I R into a singular part, 


-jkr 

>0 - S h 

e n 


n 


dA 


(3) 


and a non-singular sum of previously calculated integrals. 


*n “ *o ' 


m 

£ 

i«l 

i*n 


(4) 


The problem, then, is reduced to calculating I 


2.1 Computation of I 


Note: 


coskr n sinkr 

5 0 ■ JV ~T~ “ + 3 J’a. ~7~ dA 


e n 


e n 


( 5 ) 


but 


sin kr 

lim — 

r n -0 n 


- 4 - 

Therefore, the second integral may be evaluated numerically, if 
desired. Temporarily both will be kept together. 


For linear elements sides are straight lines: 



Figure 2. Type Q-4 Element 



Figure 3. Polar Coordinate Definition for T-3 Element 


Rephrasing I in local polar coordinates: 

e ' jlCr n B («) -Jki 

‘o ‘ ~ r n dr n dB ' J ‘o • f 0 e " dr n d ® 

e n 

where R , is the equation of line Afi in polar coordinates. 

/X Afn 





Figure 4. Polar Coordinate Definition for Q**4 Element 


For a quadrilateral element 


* Rn(®) -jfcr » R -(»*) -jkr 

L " S n S Q e dr n dS + J Q J Q 2 e dr.d®* 



Thus, an expression for R^(®) 1*2 mU8t b« * oun< * 



As shown in Figure 5, by similar triangles ABC ADE 

R 2 sin« R^coea - 
rsinO " rcosfi - R 1 

(rcosfl - Rj^HR^sina) - (rsin*) (R 2 cosa - R x ) 


r(R 2 sinacos® - R 2 cosasin« + R 1 sin») - Resina 


but, 


R 2 (sin<*cos9 - coshes in® ) - Rsin(cr -6) 


RiR^sino 




7 


or 


R l R 2 sinflt 

R " R 2 sin(a - 0) + RjSin® 

Therefore: 

R l R 2 8ina 

R l (9) “ R 2 ein(e - 0) + RjSin* 

R 2 R 3 8ini8 

R 2 (9 * ) " R 3 ein(^ - ©*) + R 2 sine* 

From Equation 6, the integral over r is non-singular and 
easily integrated. Equations 6 and 7 can now be partially 
evaluated numerically by 


r? rR(®) 
J 0 J 0 



dr n d6 



-jkr 

fi 

-j* 


n 


R<*) 

0 


60 


.if [ -jkR(e) 
k J 0 l® 


l] dfi 


A 

k 


0 


-jkR(e) 


60 


_ 11 

k 


(9) 


This integral is easily evaluated numerically using Gauss quadra- 
ture over 0. 
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3. EXTENSION TO HIGHER-ORDER ELEMENTS 

Fox many applications, linear triangular and quadrilateral 
elements are sufficient. However, for some applications, 
higher-order (e.g., quadratic or cubic) elements may provide an 
increase in accuracy. There are two types of elements with higher 
interpolation orders that are of interest here — iso-parametric 
and sub-parametric elements. 

In sub-parametric elements the geometry of the element is 
defined by a lower order of interpolation than the function 
values. If the geometry is defined in terms of linear shape 
functions the method of determining I Q is exactly the same for 
all sub-parametric elements as the method presented above, where 
m is the number of nodes in the element. 

In iso-parametric elements, however, the geometry is defined 
by higher -order interpolation, as shown in Figure 6. 



Figure 6. Isoparametric Element Displacements 


- 9 - 

Sine* the boundary of the element ia no longer dafinad by 
straight lines, tha tachniquas presented abova will not intagrata 
I correctly, as shown by tha shadad araa in tha figura. In this 
author's experience computing i for an iso-par amstric alanant is 
much mors difficult than computing it for a sub-par amstric cla- 
mant. 


In many applications tha gaomatry of tha rag ion is naithar 
complex enough to warrant iso-parametric elements nor critical 
enough that it needs to be represented exactly. In addition, 
there is some evidence which indicates that higher-order elements 
do not improve tha accuracy of tha method significantly. For 
these reasons it was decided to concentrate on tha sub-parametric 
techniques, which is much simpler and should give a good indica- 
tion of the advantages, if any, of using higher-order elements. 

If necessary, iso-parametric capability may be developed at a 
later data. 

3.1 Cpmautat i gn of the -Integral for Sub-Par amstric Elements 

At this juncture it is interesting to look at methods for 
computing I Q in general — i.e., covering the range of possible 
root nodes and element shapes. 

The root node must be one of three types: 

1. Corner node (all element types), as shown in Figure 7 

2. On a side of the element (quadratic and cubic types), 
as shown in Figure 8 
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3. In the interior of the element (certain cubic types), 
as shown in Figure 9 

In sach case the root node is marked with a dot. All other nodes 
used to define the geometry of the elements are marked with an 
"x". Letters marking each of the sub-tri angles are enclosed in 
circles. 

Note that in each case the element may be divided into from 
one to four triangles by connecting the root node with each of 
the other corner nodes in the element. The integral for the 
entire element is simply the sum of the integral for each of the 
triangles. Also note the following: 1) One or two of the trian- 

gles will collapse to lines (i.e., have area of zero) for certain 
root nodes; 2) the number of subtriangles in each element is 
equal to the number of corner nodes in the element. By using the 
test "is the area of this triangle equal to zero" it is possible 
to write a general algorithm to compute I Q for any sub-parametric 
element as shown in Table I . 
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TABLE I. 

6.1 Algorithm to Calculate 1^ 

A) Main Routine: 

1) Set the sum I to zero 

2) FOR i ■ 1 to number of corner nodes in the element (S) 
DO 

a) Select 2 corner nodes to create a triangle: 

a - i 

b - i + 1 

if (b > S ) b - 1 

b) Calculate contribution to I for this triangle 
using the routine below: 

c) Add contribution to I. 

3) END 

B) Routine to calculate contribution to I for each triangle. 

1) Calculate area of the triangle. 

2) IF (the area of the triangle is zero) THEN 

a) return a value of zero 
ELSE 

b) compute the integral over the triangle (n,a,b) 
using Eg. 9. 

3) END 



13 


4. INTEGRALS OF NORMAL DERIVATIVES 

Up to this point, only integrals of the characteristic solu- 
tion have been considered. However, it is well known that 
integrals of the normal derivative of the characteristic solution 
must also be computed. 

If an element is a portion of a plane, the normal to the 
surface is perpendicular to the vector extending from the root 
node to any point on the element. This is true for all triangles 
formed with linear sides, and is true for quadrilaterals provided 
opposite sides lie in the same plane. In such cases where the 
normal and the "position* vector are perpendicular the normal 
derivative is zero. 

Sub-parametric elements, which are being used exclusively in 
this report will be planar for all triangular elements. Sub- 
parametric quadrilateral elements may also be defined as planes 
with very little loss of versatility. Therefore, for sub- 
parametric elements, all integrals of the normal derivative are 
zero, and thus, very easy to compute. 


5. CONSTRUCTION OF A BOUNDARY ELEMENT PROGRAM 


In order to properly utilize the integration routines 
presented following this chapter, it is necessary to place them 
in the intended type of program structure. This section will 
briefly explain the type of program structure for which they were 
designed. 

The major operations of the program are: 

1. Input mesh definition into the program. 

2. Assemble the system of equations. 

3. Add internal sources to the equation system (optional). 

4. Generate additional equations for overdetermined sys- 
tem. 

5. Solve the system of equations. 

6. Extract results. 

Of course, this flow chart will vary somewhat if, for example, a 
line by line equation solver is used. 

The most complex step of the program is the assembly of the 
system of equations. This step, is discussed in more detail in 
Table 2. Figure 10 is a flow-of -control diagram, showing which 
routines are called by, or connected to other routines. 
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Figure 10 


Flow of Control for Assembling the System 
of Equations. 
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TABLE II 

Assembly of the System of Equations 

0. Define complex vectors HROW, GROW of length NOD. 

1. FOR i * 1 to the number of nodes (NOD) DO 

a. Set the vectors HROW, GROW to zero. These vectors will hold, 

.respectively, the integrals of the normal derivative of the 
characteristic solution (H) and the characteristic solution 
(G), respectively. 

b. FOR j = 1 to the number of elements (NEL) 00 

i. Localize coordinates and node numbers for element j into 
temporary arrays. 

ii. IF j is triangular THEN 

1. IF node i is in element j THEN 

call integration routine for a triangle containing 
the root node (SAINT3). f 

1, ELSE 

call integration routine for a triangle not 
containing the root node (INTE3). 

ii, ELSE 

i. If node i is in element j THEN 

call integration routine for a quadrilateral 
** containing the root node (SAINT4). 

1, ELSE 

call integration routine for a quadrilateral 
not containing the root node. 

iii. Add the integrals just computed to their proper locations 
in the vectors HROW, GROW. 

b, END of loop on the elements, j. 

c. Add the row of integrals for node i. accumulated in HROW and 
GROW to the coefficient matrix following the boundary 
conditions specified for each node. 

1, END of loop on the nodes, i. 
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The tests in l.b.ii of Table 2, which are used to select the 
type of integral to be computed, should properly be placed in a 
separate routine for easy maintenance. However, since the tests 
are quite simple and the integration routines very flexible, it 
should be more efficient to incorporate the tests directly into 
the main assembly routine. Step l.c of Table 2 is dependent upon 
the type of solution procedure chosen, and should be placed in a 
separate routine. This subject will be discussed more later. 

S.l Construction of the Integration Routines 

A listing of the integration routines (which is slightly 
out-of-date, but all that was available at the time of this writ- 
ing) is given in the Appendices, along with shape function and 
certain other utility routines. This section will briefly dis- 
cuss the structure of the integration routines. 

The routines may be divided into three groups: 

1. Computation of the I integral for an element 

COINTE, IOINTE 

2. Computation of the integrals for an entire element 
SAINT3, INTE3, SAINT4, INTE4 

3. Utility routines and shape functions SKAQUA, SHATRI, 
NORMAL along with the subroutines to SHAQUA: 

SH2DQQ, SH2DCQ 


The routines in 1 and 3 are called by the routines in 2 and are 


transparent to the remainder of the program. 


5.1.1 The Routines to Compute I Q The routines COINTE and IOINTE 
are a direct implementation of the method described previously. 
COINTE divides the sub-parametric element into the number of 
sub-triangles. IOINTE is called by COINTE to calculate the 
integral over each sub-triangle. If one angle of the sub- 
triangle is 180 degrees (which is equivalent to asking "is the 
area zero*, and is much easier to calculate), COINTE returns a 
value of zero. 

5.1.2 The Routines to Compute the Element Integrals INTE3 and 
INTE4 compute the integrals for triangles (thus the 3) and qua- 
drilaterals, respectively, for the cases when the root node is 
not in the element. These are straight forward implementations 
of standard Gaussian quadrature. 

SAINT3 and SAINT4 are images of INTE3 and INTE4 except for 
two modifications. First, during the Gaussian quadrature, the 
integral for the root node is not calculated, since it is singu- 
lar and inaccurate. Second, additional statements are added at 
the end of the routine to call COINTE to compute the I Q integral 
for the element. I Q is then used to compute the integrals for 
the root node. 

The "SAINT" and "INTE" routines are not combined because a 
large number of tests would be needed within the routines to 
decide whether the root node is in the element. These tests 
would make the computation much less efficient. 


Integration routines for triangular and quadrilateral ele- 
ments are not combined in one routine because certain operations 
may be eliminated when triangular elements are handled 
separately. Also, the dimensions of the arrays and certain DO 
loop parameters may be specified as constants, which, on some 
compilers, is more efficient than specifying variables. 

5.1.3 utility Routines The utility routines contain NORMAL, 
which calculates the normal vector for a triangle with linear 
sides, and the shape function routines. The shape functions are 
simplifications of finite element shape functions (the deriva- 
tives of the shape functions need not be calculated) for 2- 
dimensional problems. 

The shape functions have the capability to handle hierarchal 
nodes, as proposed by Zienkiewicz. Any of the midside nodes of 
the higher order elements may be removed, reducing the number of 
degrees of freedom and the order of the element. To remove node 
i, for example, specify the ith position in the element connec- 
tivity as zero, i.e., ID(i) - 0. 

The integration routines assume that the quadrilateral ele- 
ments are plane, or nearly plane. To minimize the error due to 
deviations from plane geometry, an "average” normal vector is 
used. However, results with non-plane elements may not be good. 

Prior to any call to an integration routine, the Gaussian 
quadrature points must be specified in the common blocks QPOINT 
(quadrilateral elements), TFOINT (triangular elements) and SPOINT 
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(one dimens ional Gauss points to integrate the 0 integral for 

v- 

S.2 Assembly of HROW. GROW into the Co.ff Icient Matrix 

Step l.c in Table 2 for the assembly of the system of equa- 
tions is in need of a few additional comments. As shown in the 
flow-of -control diagram, an additional routine is responsible for 
converting the integrals stored in HROW and GROW into a single 
equation for node i. The equation is then assembled into the 
system by the proper routine for the type of solution method 
used. 


The conversion of the integrals into an equation is a simple 
operation based on the boundary conditions for the nodes. The 
program BOUN3D, originally written for three dimensional heat 
transfer contains all of the instructions necessary for this 
operation. The comments in the program should be self- 
explanatory. 

Another pressing problem with boundary conditions is which 
boundary condition to enforce at a corner where different condi- 
tions are specified on each side of the corner. The system used 
in BOCJN3D, and suggested for this program is to place two (or 
more) nodes close to the corner. The different boundary condi- 
tions may be specified on different nodes. 

When the integrals are converted to an equation the correc- 
tion factor for the included angle at the node must also be 
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incorporated. It is suggested the angle be calculated from the 
geometry of the mesh, rather than the integrals , as is usually 
done for heat transfer problems. This task is probably handled 
best within a mesh generation program. Note that an additional 
array will have to be added to B0UN3D to hold the angle at each 
node. 


Finally, a note about equation solvers needs to be given. 

The usual manner of solving unsymmetric full systems by Gauss 
elimination has proven quite satisfactory in the past, and is 
suggested for this program. Iterative methods, such as SOR, BSOR 
and Gauss-Seidel usually will not converge at a reasonable rate. 

A possible variation on Gauss elimination would take advan- 
tage of the fact that boundary element equations are formed com- 
pletely, one at a time. Since the complete equation is available 
at once, it is possible to eliminate the first i-1 coefficients 
from row i at the time the row is added to the coefficient 
matrix. This "line by line" elimination would reduce the amount 
of storage needed for the matrix from N * N to N * (N -1 )/2. 
However, it would make it impossible to pivot the matrix. 

Finally, by separating the formation of the equation from 
its addition to the coefficient matrix, it will be easy to add 
specialized solution routines. For example, for some least 
squares methods, the number of equations and variables are not 
equal. Such solvers often require specialized assembly methods. 
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6. CQHCLtfglQHg 

This work has provided several insights which may be useful. 

1. Sub-parametric elements seem to be the most efficient means 
of utilizing higher-order elements. However, there is very 
little evidence to suggest that highar-order elements will 
greatly improve the accuracy of the eolution. 

2. Attention should be given to a method(s) of performing an 
iterative "design" type problem solution. It is very expen- 
sive to iterate using boundary element because the system of 
equations needs to be decomposed at each iteration. This 
could become a very pressing problem whenever iterations are 
necessary, such as a noise path identification application. 

3. For problems where all of the nodes are given the same boun- 
dary conditions (such as an impedance) there should be very 
little problem imposing boundary conditions (i.e., the prob- 
lems that haunt solid mechanics and heat transfer applica- 
tions should not appear). This is an area that can be 
ignored until a serious problem crops up 

4. The solution of the system of equations will continue to be 
the most time consuming part of the analysis. Some improve- 
ment must be made in this area, but such improvements prob- 
ably not a relevant concern at this time. 
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APPENDICES 


APPENDIX A: Subroutines COINTE and IOINTE 


function cointe (korn, iroot, x,y, z, ak) 


Computes the Io integral ovez the element by dividing the (sub-) 
parametrix element into a number of lineaz triangular elements 


Inputs: 

lroot - local number of the root node 
npe - number of nodes in the element (used to determine wheter 
element is triangular or quadrilateral) 
x,y,z - element nodal coordinates 
id - connectivity for element 
ak - constant from Helmholtx equation 


Outputs : 

cointe - value of the Io integral for the element 


Calls : 

iointe - computes Io for individual triangles in element 


Notes : 

1) iointe requires gauss integration points for one dimension. These 
must be set before routine is called. 

2) Intended only for sub-parametric elements. Serious errors may 
occuz if used with parametric elements. 


implicit real* 8 (a-h,o-r) 
complex* 16 cointe, iointe 
dimension x(korn) ,y(korn) , z(korn) 

arrays hold points for the nodes of subtriangle 

dimension xx(3) ,yy(3) ,zz(3) 

find element type, initilize variables 

cointe - (0. ,0. ) 
xx (1) - x( lroot) 
yyd) - y(iroot) 
zz(l) - z(iroot) 

divide element into subtriangles and integrate each 

il - korn 
do 100 i - l,korn 
i2 - i 

if(il .eq. lroot .or. 12 . eq. lroot) goto 100 
xx (2) - x(il) 
xx( 3) - x(i2) 
yy(2) - y( il) 
yy(3) - y(i2) 
zz(2) * z(il) 
zz(3) - z(i2) 

cointe - cointe ♦ iointe(xx,yy,zz,ak) 

00 il - 12 


J 

I 


! 

return 

end 

i 

I — — — 

function iointe(x,y,z, ak) 


Computes the constant integral over a triangular element or a triangular 
portion of an element 


Inputs: 

x,y,z - the coordinates of the triangle. The root node is the first in 
the arrays. 

ak - constant from the helmholtz equation 


Outputs : 

iointe - value of the integral of: 
exp(-jkr)/r dA 


Notes: 

1) all reals are double precision 

2) the integral is performed using both analytic and numerical techniques 

3) if the area of the triangle is zero, iointe - 0 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION X(3) ,Y(3) ,Z(3) 
complex* 16 iointe 

GAUSS POINTS POR NUMEkiCAL INTEGRATION 

COMMON /SPOINT/ ETA (16) ,W(16) ,NINPT 

TOP IS THE DOT PRODUCT OF Rl AND R2 

TOP - (X(2) - X(l)) * (X(3) - X(l) ) 

* + (Y(2) - Y(l)) * ( Y(3 ) - Y(l) ) 

* + (Z(2) - Z(l) ) * (Z(3 ) - Z(l) ) 

Rl - DSQRT( (X(2)-X(1))**2+(Y(2)-Y(1))**2+(Z(2)-Z(1))**2) 

R2 - DSQRT( (X(3)-X(1))**2+(Y(3)-Y(1))**2+(Z(3)-Z(1))**2) 

DOT R i & R2 TO FIND ANGLE OF ELEMENT 

COSA - TOP/ (Rl* R2) 

if cosa - l t the triangle is collapsed to a line 

if (dabs(dabs(cosa)-l.dO) . le. l.d-4) then 
iointe • (0.d0 r 0.d0) 
return 
end if 

compute the Io integral for the triangle 

ALPHA - DACOS (COSA) 

SINA - DS IN (ALPHA) 

NUMERATOR OF R-CARAT 

TOP - Rl * R2 * SINA 
A2 - ALPHA * 0.5D0 
iointe * (0.d0,0.d0) 

INTEGRATE THE THETA INTEGRAL 

do 100 i»l,ninpt 
RC IS r-carat 

rc - top/ (R2*dsin(a2* ( l.d0-eta( i) ) ) +Rl*dsin(a2* ( l.d0+eta( i) ) ) ) 
8 )ci , ckr are sin, cos of k* r-carat 


I skr - dsin(ak * rc) 
j ckr - dsqrt(l.dO - skr*skr) 

i multiply by weights now to avoid conversion later 

skr - skr * w(i)/ak 
ckr - ckr * w(i)/ak 

| — * Add intermdiates to the totals, GA..GC 

)0 iointe - iointe + dcmplx(skr , ckr ) 

; CORRECT FOR JACOBIAN and constant term 

iointe - iointe * dciplx(a2 , 0 -d0) - dcmplx(0.d0, alpha/ ak) 
RETURN 
END 


APPENDIX £: Subroutines INTE4, SAINT4, INTE3, and SAINT3 


subroutine inte4 ( xp , npe , x,y,z, id,g ,h,ak) 


PERFORMS INTEGRALS FOR quadrilateral elements not 
CONTAINING THE root node 


Inputs : 

XP - ARRAY OF LENGTH 3, CONTAINING THE COORDINATES OF THE NODE 
IN QUESTION (XP(I) - X(I)TH COORDINATE OF THE NODE) 
npe - number of Nodes Per Element 
x,y,z - nodal coordinates of element 

id * element connectivity (used in some shape function routines 
ak - constant from helmholtz equation 
iroot - number of the node in the element which is the root node. 
If the root is outside the element, iroot - 0 

Outputs : 

g - complex, d.p. array of length npe containing integrals 
OF U* OVER THE ELEMENT 

H - complex, D.P. array of length NPE containing integrals 
OF Q* OVER THE ELEMENT 


NOTES: 

1. THIS ROUTINE IS CAPABLE OF HANDLING any type of quadrilateral 
element for which a chape function routine is installed 

2. THE INTEGRATION POINTS GIVEN IN /POINT/ ARE FOR INTEGRATION 
OVER line. The integration points must be set before the 
routine is called. Integration is done in both coordinate 
directions in the usual manner. 

3. Maximum number of integration points presently is 16X16 

4. This routine is valid only for sub-parametric elements 

IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION XP(3),X(4),y(4),z(4), id (npe) 
complex* 16 g (npe ) , h (npe ) , us tar , qstar 

block of 'space* for all of the shape function routines 

common /sspace/ sha(12) ,xx(3) ,yy(3) ,zz(3) ,u(3) ,xc(3) 

space for integration points (1-D, gauss type) 

COMMON /QPOINT/ psi ( 16 ) , W( 16 ) , NINPT 

initialize all integrals to zero 

do 10 i - l,npe 

g ( i) - (O.dO, O.dO) 
h(i) - (O.dO, O.dO) 

Find the normal vector and area for the element 

First triangle (nodes 1, 2 & 3) 

call normal (x,y,z,u,al) 

second triangle 

xx ( 1) - x(2) 
xx(2) - x(3) 
xx(3 ) - x(*> 
yyd) - y(2) 


I! yy(2) - y(3) 

M yy(3) - y(4) 

22 ( 1 ) - 2 ( 2 ) 

22 ( 2 ) « 2 ( 3 ) 

22 ( 3 ) - 2 ( 4 ) 

call normal(xx,yy, 22 , u, a2) 
area - al + a2 

INTEGRATE! 

DO 100 I - 1, ninpt 
chi - psi(I) 
cl - l.dO - chi 
c2 - l.dO + chi 
f do 100 j - 1, ninpt 

eta - psi( j ) 
wate * w( i) * w( j ) 
el * l.dO - eta 
e2 « l.dO + eta 

compute geometric SHAPE FUNCTIONS 

sha(l) - 0. 25d0 * cl * el 

sha(2) - 0. 2SdO * c2 * el 

s ha ( 3 ) - 0. 25d0 * c2 * e2 

sha(4) - 0. 25d0 * cl * e2 

compute vector from root to integration point, (xc) 

xc(l) - - xp(l) 
xc(2) - - xp( 2) 
xc(3) - - xp(3) 
do 50 k • 1,4 

xc ( 1) - xc ( 1) + sha(k) * x(k) 

xc (2) - xc (2) + sha(k) * y(k) 

) xc(3) - xc(3) + sha(k) * 2 (k) 

R - dsqrt(xc( 1) **2 + xc(2)**2 + xc(3)**2) 

compute interpolation shape functions 

if(npe . ne. 4) call shaqua(chi, eta,npe, id, sha) 

dot normal, u, and vector, r, to get angle 

cosa - (xc(l)*u(l) + xc(2)*u(2) + xc(3)*u(3) )/r 
ustar - dcmplx(dcos(ak * r ) *wate/r , -dsin(ak*r ) *wate/r ) 
qstar » ustar * dcmplx(-cosa/r ,-ak*cosa) 

DO ADDITIONS FOR INTEGRATIONS 

DO 100 k * l,npe 

H(k) - H(k) + dcmplx(SHA(k) , O.dO) * QSTAR 
G(k) - G(k) + dcmplx(SHA(k) , O.dO) * USTAR 
)0 continue 

CORRECT FOR JACOBIAN (make temporary use of ustar) 

ustar - dcmplx(area*0. 25d0 , O.dO) 

DO 110 I - l,npe 

0(1) - G(I) * ustar 
L0 H(I) • H ( I ) * ustar 

RETURN 
END 
$ EJECT 

SUBROUTINE saint4 (npe , x,y , z , id,g,h,ak, iroot) 


j PERFORMS INTEGRALS FOR quadrilateral elements 
| CONTAINING THE root NODE 

I 

j Inputs: 

f npe - number of Nodes Per Element 
( x,y,z - nodal coordinates of element 

id - element connectivity (used in some shape function routines 
\ ak - constant from helmholtz equation 

[ iroot - number of the node in the element which is the root node. 
If the root is outside the element, iroot * 0 


i Outputs: 

g - complex, d.p. array of length npe containing integrals 
OF U* OVER THE ELEMENT 

H - complex, D.P. array of length NPE containing integrals 
OF Q* OVER THE ELEMENT 


NOTES: 

1. THIS ROUTINE IS CAPABLE OF HANDLING any type of subparametr ic 
quadrilateral elements which a shape function routine is 
installed 

2. THE INTEGRATION POINTS GIVEN IN /POINT/ ARE FOR INTEGRATION 
on a line. These are used in both directions for integrating 
over the element. (Must be set before routine is called) 

3. This routine may call function iointe, which needs 1-D gauss 
points in /spoint/. These must be set before routine is called 

4. Only valid for planar, subparametr ic elements 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION XP(3),X(4) ,y(4) ,z(4) , id(npe) 
complex* 16 g(npe) ,h(npe) , ustar , cointe 

block of 'space* for all of the shape function routines 

common /sspace/ sha(12) ,xx(3) ,yy(3) ,zz(3) ,u(3) ,xc(3) 

space for integration points (1-D, gauss type) 

COMMON /QPOINT/ psi(16) ,W(16) ,NINPT 

initialize all integrals to zero 

do 10 i * l,npe 

g(i) « (0. dO , 0 .d0) 
t h(i) - (0.d0,0.d0) 

Area for the element: for first triangle 

call normal(x,y,z,u,al) 

Becond triangle 

xx (1) - x(2) 
xx (2) - x( 3) 
xx (3) - x(4) 
yy(i) - y(2) 
yy(2) - y(3) 
yy(3) - y(4) 
zz(l) - z(2) 
zz(2) - z(3) 
zz(3) - z(4) 

call normal (xx,yy,zz,u,a2) 

area - al + a2 


f INTEGRATE! 

DO 100 I - l,ninpt 
chi - psi(I) 
cl - l.dO - chi 
c 2 - l.dO + chi 
j do 100 j • 1 , ninpt 
eta - psi(j) 
wate - w( i) * w(j ) 
el » l.dO - eta 
; e2 » l.dO + eta 

!■ compute geometric SHAPE FUNCTIONS 

sha(l) - 0. 25d0 * cl * el 

sha(2) - 0. 25d0 * c2 * el 

aha(3) - 0.25d0 * c2 * e2 

sha(4) - 0. 25d0 * cl * e2 

compute vector from root to integration point , (xc) 

xc(l) - - xp( 1) 

xc(2) - - xp(2) 

xc(3) - - xp(3) 

do 50 k - 1,4 

xc(.l) - xc(l) + sha(k) * x(k) 

xc(2) - xc(2) + sha(k) * y(k) 

xc(3) - xc(3) + sha(k) * z(k) 

R - dsqrt(xc( 1) **2 + xc(2)**2 + xc(3)**2) 

us tar - dcmplx(dcos (ak * r)*wate/r,-dsin(ak*r)*wate/r) 

Compute interpolation shape functions 

if(npe .ne. 4) call shaqua(chi,eta,npe, id,sha) 

DO ADDITIONS FOR INTEGRATIONS 

DO 100 JJ - l,npe 

skip the root node, if one exists 

if(jj .ne. iroot) then 

G(JJ) - G(JJ) + dcmplx(SHA( JJ) ,0.d0) * USTAR 
end if 
i0 continue 

CORRECT FOR JACOBIAN (make temporary use of ustar) 

ustar * dcmplx(area * 0.25d0,0.d0) 

DO 110 I - l,npe 
.0 G(I) - G(I) * ustar 

if root is in the element correct for the constant term 

if (iroot .ne. 0) then 
korn - 4 

g( iroot) - cointe(korn, iroot, x,y,z,ak) 
do 200 i - l,npe 

if(i .ne. iroot) g( iroot) ■ g( iroot) - g(i) 

10 continue 

end if 

RETURN 

END 
1 EJECT 

SUBROUTINE INTE3(XP,npe,X,y,z, id,G,H,ak) 


PERFORMS INTEGRALS FOR subparametr ic triangular elements not 


CONTAINING THE root node 


Inputs: 

XP - ARRAY OP LENGTH 3, CONTAINING THE COORDINATES OP THE NODE 
IN QUESTION (XP(I) - X(I)TH COORDINATE OP THE NODE) 
npe - number of Nodes Per Element 
x,y,z - nodsl coordinates of element 

id - element connectivity (used in some shape function routines 
ak - constant from helmholtz equation 


! Outputs : 

g - complex, d.p. array of length npe containing integrals 
OF U* OVER THE ELEMENT 

H - complex, D.P. array of length NPE containing integrals 
OF Q* OVER THE ELEMENT 


NOTES: 

1. THIS ROUTINE IS CAPABLE OF HANDLING any type of subparametr ic 
triangular element 

2. THE INTEGRATION POINTS GIVEN IN /POINT/ ARE POR INTEGRATION 
for a triangular region. The integration points must be set 
before the routine is called. 

3. Maximum number of integration points presently is 16 

4. This routine is valid only for sub-par ametr ic elements 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION XP(3) ,X(3) ,y(3) ,z(3) , id(npe) 
complex* 16 g(npe) ,h(npe) , us tar ,qstar 

• — block of 'space' for all of the shape function routines 

common /sspace/ sha(12) ,xx(3) ,yy(3) ,zz(3) ,u(3) ,xc(3) 

space for inte cation points (i-D, gauss type) 

COMMON /TPOINT/ all(16) ,al2(16) ,W(16) ,NINPT 

initialize all integrals to zero 

do 10 i - l,npe 

g(i) - (O.dO, O.dO) 
h(i) - (O.dO, O.dO) 

Find the normal vector and area for the element 

call normal(x,y,z,u,area) 

INTEGRATE! 

DO 100 I - l,ninpt 

use shape functions to caluculate vector from root 

to point in the element, (xc) 

xc ( 1) - all (i)*(x(l)-x(3)) + al2 (i)*(x(2)-x(3)) + x(3) - xp(l) 

xc(2) - all(i)* (y(l)-y(3) ) + al2(i)*(y(2)-y(3)) + y(3) - xp(2) 

xc(3) - all(i)*(z(l)-z(3)) + al2(i)*(z(2)-z(3)) + z(3) - xp(3) 

R - dsqrt(xc(l) **2 + xc(2)**2 + xc(3)**2) 

compute interpolation shape functions 

call sh&tr i ( all ( i ) , al2 ( i ) , npe , id , sha ) 

dot normal, u, and vector, r, to get angle 

cosa - (xc(l)*u(l) + xc(2)*u(2) + xc(3) *u(3) )/r 
ustar - dcmplx(dcos(ak * r)*w(i)/r,-dsin(ak*r)*w( i)/r) 
qstar * ustar * dcmplx ( -cosa/r ,-ak*cosa) 

DO ADDITIONS POR INTEGRATIONS 


DO 100 k - l,npe 

H|,K) - H(k) + dcmplx(SHA(k) , O.dO) * QSTAR 
G(k) - G(k) + dcmplx(SHA(k) , O.dO) * USTAR 
0 continue 

CORRECT FOR JACOBIAN (make temporary use of ustar) 

ustar - dcmplx( area, O.dO) 

DO 110 I - l,npe 

G( I ) - G ( I ) * ustar 
.0 H(I) - H(I) * ustar 

RETURN 

END 

SUBROUTINE saint3 (npe , x ,y , z , id,g,h,ak, iroot) 


PERFORMS INTEGRALS FOR aubparametr ic triangular elements not 
CONTAINING THE root node 


Inputs: 

npe - number of Nodes Per Element 
x,y,z - nodal coordinates of element 

id - element connectivity (used in some shape function routines 
ak - constant from helmholtz equation 
iroot - local number of the root node 


Outputs : 

g - complex, d.p. array of length npe containing integrals 
OF U* OVER THE ELEMENT 

H - complex, D.P. array of length NPE containing integrals 
OF Q* OVER THE ELEMENT 


NOTES: 

1. THIS ROUTINE IS CAPABLE OF HANDLING any type of aubparametr ic 
triangular element 

2. THE INTEGRATION POINTS GIVEN IN /POINT/ ARE FOR INTEGRATION 
for a triangular region. The integration points must be set 
before the routine is called. 

3. Maximum number of integration pointB presently is 16 

4. This routine is valid only for sub-* parametric elements 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION XP(3),X(3),y(3),z(3), id (npe) 
complex* 16 g (npe) ,h(npe) , ustar , cointe 

block of 'space' for all of the shape function routines 

common /sspace/ sha(12) ,xx(3) ,yy(3) , zz(3) ,u(3) ,xc(3) 

space for integration points (1-D, gauss type) 

COMMON /TPOINT/ all ( 16 ) , al2 ( 16 ) , W( 16 ) , NINPT 

initialize all integrals to zero 

do 10 i - l,npe 

g ( i ) - (O.dO, O.dO) 
h(i) - (O.dO, O.dO) 

Find the normal vector and area for the element 

call normal (x,y,z,u, area) 

INTEGRATE! 


D 




DO 100 I « l,ninpt 

t use shape functions to caluculate vector from root 

to point in the element, (xc) 

xc(l) - all(i)*(x(l)-x(3) ) + al2 (i)*(x(2)-x(3)) + x(3) - xp(l) 

xc(2) - all( i)*(y(l)-y(3) ) ♦ al2( i)*(y(2)-y(3) ) + y(3) - xp(2) 

xc(3) - all( i) * (z( l)-z(3) ) + al2(i)*(z(2)-z(3)) + z(3) - xp(3) 

R - dsqrt(xc(l)**2 + xc(2)**2 + xc(3)**2) 

compute interpolation shape functions 

call shatr i ( all ( i ) ,al2(i) ,npe, id,sha) 

dot normal, u, and vector, r, to get angle 

ustar - dcmplx(dcos(ak * r ) *w( i)/r , -dsin(ak*r ) * 1 w( i)/r ) 

DO ADDITIONS FOR INTEGRATIONS 

DO 100 k * l,npe 

if (k .ne. iroot) G(k) • G(k) + dcmplx(SHA(k) , O.dO) * USTAR 
0 continue 

CORRECT FOR JACOBIAN (make temporary use of ustar) 

ustar • dcmplx( area, O.dO) 

DO 110 I » l,npe 
0 G(I) - G(I) * ustar 

correct for root node 

g( iroot) - cointe(npe,x,y,z,ak) 
do 120 i - l,npe 

0 if(i .ne. iroot) g( iroot) - g( iroot) - g(i) 


RETURN 

END 


PENDIX C: Subroutines NORMAL, SHAQUA SHZDQQ, SHATRI , and SH2DCQ 


subroutine normal (x,y,z,u, area) 


Computes the normal and area of a triangular plane in three 
dimensions 


Inputs: 

x,y,z - coordinates of the three nodes 

u - normalized normal vector formed by crossing the vector 
between nodes 1 and 2 into the vector between nodes 
1 and 3 

area - area of the triangle 

COMPUTE VECTORS R1 £ R2 ALONG SIDE 1 £ 3 OF ELEMENT 

implicit real*8 (a-h,o-z) 

dimension x(3),y(3),z(3),u(3),rl(3),r2(3) 

form the two vectors to be crossed 

rl(l) - x(2) - x ( 1 ) 

rl(2) - y ( 2 ) - y(l) 

rl(3) - z ( 2) - z ( 1) 

r2(l) - x(3) - x(l) 

r2(2) - y(3) - y(l) 

r2(3) - z(3) - z ( 1) 

CROSS R1 & R2 TO FIND U, NORMAL TO THE SURFACE 

U(l) - (Rl( 2) * R2(3) - Rl( 3 ) * R2(2)) 

U( 2) - (Rl(3) * R2(l) - Rl( 1) * R2 ( 3 ) ) 

U(3) - (Rl(l) * R2(2) - Rl( 2) * R2(l)) 

THE MAGNITUDE OF U IS TWICE THE AREA OF ELEMENT 

BOT - DSQRT (U(l)**2 + U(2)**2 + U(3)**2) 

AREA • 0.5 * BOT 

NORMALIZE U 

U(l) - U( 1) /BOT 
U(2) - U(2)/BOT 
U(3) - U ( 3 ) /BOT 

return 

end 


MODULE FOR 2-D planar shape functions 


SUBROUTINE SHagua(SS ,TT,npe, id,SHA) 


CONTROL FOR TWO-DIMENSIONAL ISOPARMETRIC ELEMENTS 
TYPES SUPPORTED: 

NUMBER HIERARCH I AL TYPE OF ELEMENT 

OF NODES NODES POSSIBLE 


4 

8 YES 


LINEAR QUADRILATERAL 
QUADRATIC QUADRILATERAL 


YES 

YES 


LAGRANGE QUADRILATERAL 
CUBIC QUADRILATERAL 


; 9 
| 12 


INPUT ARGUMENTS: 

SS,TT — NATURAL COORDINATES 

npe — number of nodes in element 
ID — element connectivity 


Outputs : 

SHA -- SHAPE FUNCTIONS AT (SS,TT) 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION S(4) ,T(4) , SHA (npe) 

DATA S/-0 . 5D0 , 0 . 5D0 , 0 . 5D0 , -0 . 5D0/ , T/-0 . 5D0 f -0 . 5D0 f 0 . 5D0 , 0 . 5D0/ 

FORM 4-NODE QUADRILATERAL SHAPE FUNCTIONS 

DO 100 1-1,4 

D SHA( I )- ( 0 . 5+S ( I ) *SS ) * ( 0 . 5 + T(I)*TT) 

ADD higher order terms if necessary 

IF (NPE .EQ. 8 .or. npe .eq. 9) CALL SH2DQQ(SS,TT,SKA, ID, NPE) 
IF (NPE .EQ. 12) CALL SH2DCQ(SS ,TT,SHA, ID, NPE) 

RETURN 

END 

EJECT 


SUBROUTINE SH2DQQ(S,T,SHA, IX, NPE) 
IMPLICIT REAL* 8 (A-H,0-Z) 
DIMENSION IX(npe) , SHA (npe) 

S2 - 0. 5d0 * (l.dO - S*S) 

T2- 0. 5d0 * ( l.dO - T*T) 

DO 100 I -5, npe 

0 SHA(I) - 0. 


if ( ix( 5) 

.ne. 

0) 

SHA( 5) - 

S2 

* 

(i.< 

if ( ix(6) 

.ne. 

0) 

SHA( 6 ) - 

T2 

* 

d-< 

if ( ix(7) 

.ne. 

0) 

SHA( 7 ) - 

S2 

* 

d.< 

if (ix(8) 

.ne. 

0) 

SHA( 8 ) - 

T2 

* 

(i.< 

if (npe . 

eq. 9) 

then 




if (ix(9) 

.ne. 

0) 

then 




SHA( 9) 

- 4. 

* 

S2 * T2 





MIDSIDE NODES 
10 - T) 

10 + S) 

10 + T) 

10 - S) 

LAGRANGE HEREAFTER 


CORRECT EDGE FOR INTERIOR 


DO 105 1-1,4 

SHA(I) - SHA(I) - 0.25 * SHA(9) 

5 IF ( IX( 1+4) .NE. 0) SHA( 1+4) - SHA(I+4) - .5 * SHA(9) 

end if 


end if 
7 K-8 


CORRECT CORNER FOR MIDSIDE 


DO 109 1-1,4 
L-I+4 

SHA(I) - SHA(I) - . 5* (SHA(K) <1- SHA(L) ) 





9 


K-L 
RETURN 
END 
EJECT 

SUBROUTINE SHATRI (Al,A2,NPE, IX, SHA) 

IMPLICIT REAL* 8 (A-H,0-Z) 

REAL* 8 L(3) 

DIMENSION SHA(npe) ,S2(3) ,S3 (3) , IX(npe) 

DATA S2/1.D0 , 0 .DO , -1.D0/ ,S3/0 .DO , 1.D0, -1.D0/ 

FORM THIRD SHAPE FUNCTION 

L(l) - A1 
L(2) - A2 

L(3) - 1. - L( 1) - L(2) 

FORM LINEAR SHAPE FUNCTIONS 

DO 10 1-1,3 

SHA( I ) - L( I ) 

IF (NPE .NE. 6) RETURN 

DO 20 J-4,6 

0 SHA(J) - O.dO 

FORM QUADRATIC TERMS AS NECESSARY 

if ( ix(4) .ne. 0) SHA(4) - 4.d0 * L(l) * L(2) 

if ( ix( 5) .ne. 0) SHA(5) - 4.d0 * L(2) * L(3) 

if ( ix( 6) .ne. 0) SHA(6) - 4.d0 * L(3) * L(l) 

— correct corners for midside nodes 

KK - 6 
DO 60 1-1,3 
JC-I+3 

SHA(I) - SHA(I) - 0 . 5d0* (SHA(KK) + SHA(K) ) 

1 KK-K 

RETURN 

END 

EJECT 

SUBROUTINE SH2DCQ(S,T,SHA, IX, NPE) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION IX(npe) , SHA(npe) 

1 2 3 4 5 

DATA S 1/0 . DO , 0 . DO , 0 . DO , 0 . DO , - 1 . DO 

6 7 8 9 10 11 12 

* , 1.D0 , 1.D0 , 1.D0 , 1.D0 , -1.D0 , -1.D0 , -1 .DO/ 

1 2 3 4 5 6 

DATA Tl/0 . DO , 0 . DO , 0 . DO , 0 . DO , ~1 . DO , ~1 . DO , 

7 8 9 10 11 12 

* -1 . DO , 1 . DO , 1 . DO , 1 . DO , 1 . DO , -1 . DO/ 


c 

- 

9. dO/32. dO 

ss 

- 

l.dO - S*S 

TT 

m 

l.dO - T*T 

S2 

m 

0 . 5d0 * SS 

T2 

m 

0. 5d0 * TT 

DO 

10 

i 1-5,12 



3HA(I)-0.d0 


ALLOW FOR MISSING NODES 


i 


IF ( I X ( S ) .ne. 0) then 
IF ( IX(6) . ne. 0) then 

SHA( 5 ) - C*(1.<S0 - T) * SS * (l.dO 
SHA(6) - C*(l.dO - T) * SS * (l.dO 
SKA( 1) - SKA(l) - (2.d0 * SHA( 5) ♦ 
SHA( 2 ) - SHA( 2 ) - (2.dO * SHA(6) + 

else 

SHA( 5 ) - S2 * (l.dO - T) 

SHA(l) - SKA(l) - O.SdO * SHA( 5) 
$HA(2) - SHA( 2) - 0. 5d0 * SHA(5) 
end if 
end if 

— side 2, nodes 7 end 8 

IF ( IX( 7 ) .ne. 0) then 
IF ( I X ( 8 ) .ne. 0) then 

SHA(7) - C*(l.dO ♦ S) * TT * (l.dO 
SHA(8) - C* (l.dO + S) * TT * (l.dO 
SHA(2) - SHA( 2) - (2.d0 * SHA(7) + 
SHA( 3 ) - SHA(3) - (2.d0 * SHA(8) + 

else 

SHA(7) - T2 * (l.dO ♦ S) 

SHA(2) - SHA( 2) - U.SdQ * SHA(7) 
SHA( 3 ) - SHA( 3) - 0.5d0 * SHA(7) 
end if 
end if 

side 3, nodes 9 end 10 

IF ( IX( 9 ) .ne. 0) then 
IF ( IX( 10) .ne. 0) then 

SHA(9) - C* (l.dO ♦ T) * S3 * (l.dO 
SHA(IO) - C* (l.dO + T) * SS * (l.dO 
SHA( 3 ) - SHA(3) - (2.d0 * SHA(9 ) ♦ 
SHA(4) - SHA(4) - (2.d0 * SHA(IO) + 

else 

SHA(9) - S2 * (l.dO ♦ T) 

3HA( 3 ) - SHA( 3) - 0.5d0 * SHA(9) 

SHA(4) - SHA(4) - 0.5d0 * SHA(9) 

end if 
end if 

side 4, nodes 11 end 12 

IF(IX(11) .ne. 0) then 
IF ( IX( 12) .ne. 0) then 

SHA(ll) - C* ( l.dO - S) * TT * (l.dO 
SHA( 12) - C*(l.dO - S) * TT * (l.dO 
SHA(4) - SHA(4) - (2.d0 * SHA(ll) + 
SHA(l) - SHA(l) - (2.d0 * SHA(12) + 

else 

SHA(ll) - T2 * (l.dO - S) 

SKA ( 4 ) - SHA(4) - O.SdO * SHA(ll) 

SHA(l) - SHA(l) - O.SdO * SKA ( 11) 

end if 
end if 
return 
END 


- 3 .dO*S) 

+ 3 .dO*S) 
SHA(6) )/3.d0 
SHA(S) )/3.d0 


• 3.dO*T) 

+ 3.dO*T) 
SHA(8) )/3.d0 
SHA(7) )/3.d0 


♦ 3 .dO * S) 

- 3 .dO * S) 
SHA(10))/3.d0 
SHA( 9 ))/3.dO 


+ 3.d0 * T) 

- 3.d0 * T) 
SHA(12) )/3.d0 
SHA(ll))/3.dO 


